
clear all;
clc;
disp('test result with Z formula using mean and sigma until the current point time');
load dataAnalysisQ_G1000k.mat

INSAMPLE_YR = 2008;
INSAMPLE_QR = 4;

inSampleLength = (INSAMPLE_YR - 1998)* 4 + INSAMPLE_QR;
outSampleLength = length(SVD)-inSampleLength;

temp =zeros(NUMBER_OF_YEAR * Quarter_LIST_LENGTH,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);

avgRollingTM =zeros(outSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH); 
bin=zeros(outSampleLength,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH); 
PHat=zeros(RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
I=ones(1,1,9,10);
filename='optimization_07q4_jun2_2.xls';

% stdProbQ
% stdProbAvgQ
% defaultFreq


%To calculate the credit score/bins----------------------------------------



    for i = 1: NUMBER_OF_YEAR
        for j = 1: Quarter_LIST_LENGTH
            temp((i-1)*4 + j, :, :) = durationCount(i,j,:,:);
        end
    end


for t = 1: outSampleLength
    xx = zeros (RATE_LIST_LENGTH-1,RATE_LIST_LENGTH); 
    
for tt = 1 : inSampleLength + t -1
    for i=1:RATE_LIST_LENGTH-1
        for j=1:RATE_LIST_LENGTH
            x(:,:) =  temp(tt, :, :);
            xx = x + xx;
        end
    end
end

for i = 1: RATE_LIST_LENGTH-1
   avgRollingTM (t, i, :) = xx(i,:)./ sum(xx(i,:),2);        
end


for i=1:RATE_LIST_LENGTH-1
    for j=1:RATE_LIST_LENGTH
        cdf = sum( avgRollingTM(t, i, j:RATE_LIST_LENGTH) );
        bin(t,i,j)=norminv(cdf,0,1);
    end
end
   
end

disp('finish Rolling Avg Transition Matrix and Bin');



% for i=1:RATE_LIST_LENGTH-1
%     for j=1:RATE_LIST_LENGTH
%         
%         % without this if, the first cell of 2nd and 3rd row cannot be 'inf' 
%         if (i == 2 && j== 1) || (i == 3 && j== 1)
%             cdf = round (sum(stdProbAvgQ(i,j:RATE_LIST_LENGTH)));
%         else cdf = sum(stdProbAvgQ(i,j:RATE_LIST_LENGTH));
%         end
%         
%         bin(t,i,j)=norminv(cdf,0,1);
%     end
% end

% xlswrite(filename, bin,'bins');


%To calculate credit cycle index Z, real and forecast----------------------

Z = -zscore(SVD);

% Website tutoring upload data from xls to matlab:
% http://blinkdagger.com/matlab/matlab-using-xlsread-to-import-excel-data/

[fitSVD] = xlsread('fittedReg_Sample98Q1_08Q4.xlsx','B6:B49');
[fcstSVD] = xlsread('fittedReg_Sample98Q1_08Q4.xlsx','C50:C53');


Zhat = zeros (outSampleLength, 1);

for i = 1 : outSampleLength
SVDhat = SVD (1:inSampleLength + i -1, :);

SVDhat(inSampleLength + i, :) = fcstSVD (i, :);
ZhatVector = - zscore (SVDhat);
Zhat(i, :) = ZhatVector(size(ZhatVector,1),:);

% Zhat(i,:) = - ( fcstSVD(i,:)- mean(SVD(1:inSampleLength,:)) )/std(SVD(1:inSampleLength,:));
end

% To minimize the distance to get estimate of w----------------------------

W = zeros (outSampleLength, 3); % first column is for metric BFS, the second colum is for metric D2
PHat_BFS = zeros (outSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
PHat_D2 = zeros (outSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
PHat_SVD = zeros (outSampleLength, RATE_LIST_LENGTH-1,RATE_LIST_LENGTH);
R = RATE_LIST_LENGTH-1;
C = RATE_LIST_LENGTH;
yrEst= INSAMPLE_YR - 1998 + 1 + floor(INSAMPLE_QR / 4);

for i = 1: outSampleLength

qrEst= mod(INSAMPLE_QR,4)+ i;
N(:,:) = durationCount(yrEst,qrEst,:,:);
P(:,:) = stdProbQ(inSampleLength+i,:,:);
X(:,:) = bin(i,:,:);
wantedZ = Zhat(i); 

% BFS distance
obj1=@(w)difBFS(w,N,P,X,wantedZ,R,C); 
[foundW,fval,exitflag,output] = fminbnd(obj1,0,0.5);
W(i,1) = foundW;
% figure;
% subplot(outSampleLength,1,i); fplot(obj1,[-0.5 0.5],'*')
% title('Mininization Object with BFS distance',... 
%   'FontWeight','bold')
PHat_BFS(i,:,:) = calculatePHat(foundW,X,wantedZ,R,C);
% sheetlist=strcat('',year_list_text(yrEst,:),'.',quarter_list_text(qrEst,:),'');
% xlswrite(filename, PHat, sheetlist);

% D2 distance
obj2=@(w)difD2(w,P,X,wantedZ,R,C); 
[foundW,fval,exitflag,output] = fminbnd(obj2,0,0.5);
W(i,2) = foundW;
PHat_D2(i,:,:) = calculatePHat(foundW,X,wantedZ,R,C);


% SVD
obj3=@(w)difSVD(w,P,X,wantedZ,R,C); 
[foundW,fval,exitflag,output] = fminbnd(obj3,0,0.5);
W(i,3) = foundW;
PHat_SVD(i,:,:) = calculatePHat(foundW,X,wantedZ,R,C);
% xx(:,:) = PHat_SVD(i,:,:);
% sheetlist=strcat('',year_list_text(yrEst,:),'.',quarter_list_text(qrEst,:),'');
% xlswrite(filename, xx, sheetlist);
% clear xx;
end




% MAE ---------------------------------------------------------------------


% naive1=stdProbAvgQ;
% naive2=stdProbQ(inSampleLength+i,:,:);
% % n(:,:)=N(yrp,freqjp,2:10,13);
% % nn=repmat(n,1,10);
% % nt(1,1,:,:)=nn;
% % Preal=P(:,:,2:10,2:11);
% e_mod=(nt.*(Preal-PZ).^2)./(PZ.*(I-PZ));
% e_avg=(nt.*(Preal-naive1).^2)./(naive1.*(I-naive1));
% e_pre=(nt.*(Preal-naive2).^2)./(naive2.*(I-naive2));
% perf_mod = mae(e_mod);
% perf_naive1 = mae(e_avg);
% perf_naive2 = mae(e_pre);
% 
% 
% Parameters={'foundW','Zhat', 'perf_mod', 'perf_naive1','perf_naive2'; foundW Zhat perf_mod perf_naive1 perf_naive2}';
% xlswrite(filename, Parameters,'Parameters');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Naive1 MAE of L1,L2 NSD,D1,D2   %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% L1_n1=Preal-naive1;
% L2_n1=(Preal-naive1).^2;
% NSD_n1=((Preal-naive1).^2)./(naive1);
% % D1_n1=
% % D2_n1=
% 
% perf_naive1e_L1=mae(L1_n1);
% perf_naive1e_L2=sqrt(mae(L2_n1).*90)/90;
% perf_naive1e_NSD=mae(NSD_n1);
% % perf_naive1e_D1=
% % perf_naive1e_D2=
% 
% disp('fini');

 
